#ifndef LYXUTILS_H
#define LYXUTILS_H
//explanation:some common utility functions
//在VS中设置：项目-属性-配置属性-C/C++ -预处理器-预处理器定义中要添加_CRT_SECURE_NO_WARNINGS,否则lyxutils::debug::getTimeString函数会出现线程不安全报错
#include <string>
#include <iostream>
#include <ctime>
#include <fstream>
#include <sstream>
#include <vector>
//#include <queue>
#include <list>
#include <map>
#include <iomanip>
#include <algorithm>
#include <cmath>
#include <thread>
#include <functional>

// uncomment following line to enable deprecated namespace algorithm(not suggested)
//#define enable_deprecated

/**\ self made lyxutils package. For convenience recommend giving an alias to the lengthy namespaces.
 * \also: aliases of the namespaces help to reduce the impact of any possible
 * \        future changes of names of the namespaces defined here
 * \usage: #include <Utils.h>
 * \       namespace dbg = lyxutils::debug;     //utilities for debug
 * \       namespace str = lyxutils::str_utils; //utilities for string operations
 * \       namespace xio = lyxutils::io;        //utilities for io operation(including file io and std io)
 * \       namespace np  = lyxutils::numpy;     //utilities for numpy-like operations(not vectorized)
 * \       namespace xos = lyxutils::os;        //utilities for os operation(ls file,change directory and path operations,etc)
 * \       namespace geo = lyxutils::geometry;  //utilities for geometry data structure and algorithms
 * \       namespace color=lyxutils::colors;    //color utilities
 * \note: if don't want to define aliases for the above namespaces but want to use the abbreviations,just:
 *         using namespace lyxutils::abbreviations;   //to use default abbreviations for all namespaces
 * \note: functions in namespaces start with double underscore(`__`) are not supposed to be used
  */
namespace lyxutils
{
    /**\brief encapsulate some functionalities for debug
     */
	namespace debug
	{
		/**\brief get the current time string
		  *\return the current time string in yyyy-mm-dd HH:MM:SS format
		  */
		std::string getTimeString();
		/**\brief time the function f
		  *\param[in] f the function to be timed
		  *\note: the unit of time is millisecond,the function being timed should take no argument
		  *\input parameter or should be adapted using lambda function or bind function
		  *\function can be easily adapted using lambda function, for example,function like:
		  *\int f(params) can be adapted like time=timeIt([&result,params](){result=f(params);});
		  */
		template<class Func>
		double timeIt(Func f)
		{
			double start, stop, duration;
			start = clock();
			f();
			stop = clock();
			duration = stop - start;
			return duration;
		}
		/**
		 * make the current thread sleep for millisecs(in milliseconds)
		 * @param millisecs time duration to sleep for
		 */
		void sleep(int millisecs);
		/**
		 * estimate total rows in a given file. if the file open failed,throw runtime_error
		 * @param file file name
		 * @param pre_read number of pre-read lines
		 * @return total estimated lines in the file
		 */
		size_t estimate_file_rows(std::string file, int pre_read=30);
	}
	/**\brief encapsule some convenient string operation
	  */
	namespace str_utils
	{
        const int __MIN_WORDWRAP_WIDTH=5;
		/**\brief split the string with given delimiter
		*\param[in] str the string to be splitted
		*\param[in] delimiter the delimiter(note:can only be single char, if several are given, the delimiter is any one of them)
		*/
		std::vector<std::string> split(const std::string &str, const std::string &delimiter);
		/**\brief split the string with a string
		  */
		std::vector<std::string> split_with_string(const std::string &str, const std::string &delString);
		/**
		 * split the string with a string,this function will always return n+1 sub strings if str contains n delString
		 * this version will retain empty string,eg "google" split with "o" will be ["g","","gle"]
		 * while split_with_string version will generate result:["g","gle"]
		 * @param str
		 * @param delString
		 * @return
		 */
		std::vector<std::string> split_with_string_re(const std::string &str, const std::string &delString);
		/**\brief to lower case
		  */
		std::string lower(const std::string &str);
		/**\brief to upper case
		  */
		std::string upper(const std::string &str);
		/**\brief replace old string with new string
		 * \param[in] str original string
		 * \param[in] oldStr,newStr old string to be replaced and new string to replace
		 * \param[in] count replace first count oldStr(default:-1 means no limit)
		  */
        std::string replace(const std::string &str, const std::string &oldStr, const std::string &newStr, int count=-1);
        /**\brief replace oldStr to newStr in str from position start to end
         *
         * @param str original string
         * @param oldStr string to be replaced
         * @param newStr string to replace
         * @param start start position, negative index means count from tail to head
         * @param end end position, negative index means count from tail to head
         * @return new string with oldStr replaced by newStr
         */
        std::string replace(const std::string &str, const std::string &oldStr, const std::string &newStr, int start, int end);
        /**\brief convert type value to string, type should implement << operator
         */
        template<typename T>
        std::string type2str(const T &t){
            std::stringstream ss;
            ss<<t;
            return ss.str();
        }
        /**\brief convert string to type value, type should implement >> operator
         */
        template<typename T>
        T str2type(const std::string &str){
            std::stringstream ss(str);
            T t;
            ss>>t;
            return t;
        }
        /**\brief concatenate element from begin iterator to end iterator using link
         */
        template<typename iterator>
        std::string join(const std::string &link,iterator begin,iterator end){
            std::string result;
            for(auto it=begin;it!=end;++it){
                if(it!=begin)result+=link;
                result+=type2str(*it);
            }
            return result;
        }
        /**\brief count occurance of sub in str from index start to end
         */
        int count(const std::string &str,const std::string &sub,int start,int end);
        /**\brief repeat str count times
         */
        std::string multiply(const std::string &str,int count);
        /**\brief center the str and fill space at both ends with c
         */
        std::string center(const std::string &str, int width, char c);
        /**\brief left align the str and fill space at right end with c
         */
        std::string left(const std::string &str,int width,char c);
        /**\brief right align the str and fill space at left end with c
         */
        std::string right(const std::string &str,int width,char c);
        /**\brief auto wrap word
         */
        std::vector<std::string> word_wrap(const std::string &str,int width);
        /**\brief frame the text like bellow
         *
         *          **************This is a title******************
         *          *                                                                                        *
         *          *         here goes some explanation or p-          *
         *          *         atent.This function deals with                  *
         *          *         word wrap and alignment,etc.wi-           *
         *          *         -th this method it's easy to g-                    *
         *          *         -enerate a framed text in a pr-                   *
         *          *         -ogram.                                                                *
         *          *                                                                                         *
         *          *********************************************
         *|------|------|--------------------------------------------|------|
         * shift  pad                          text area                            pad
         *|------|----------------------------------------------------------|
         *                                                 width
         */
        std::string frame(const std::string &title,const std::string &text,int width,int pad,char border,int shift=0);
        /**\brief strip the chars at both ends
         */
        std::string strip(const std::string &str,const std::string &chars=" ");
		/**\brief find the first occurrence of sub in str. if not found return -1
		 */
		int find(const std::string &str,const std::string &sub,int start=0);
//        template<typename T,typename allocator=std::allocator<T> >
//        void literal_range(std::string literal,std::vector<T,allocator> &range){
//            std::vector<std::string> paras=split_with_string_re(literal,":");
//            if(paras.size()==1){
//                range.push_back(str2type<T>(paras[0]));
//                return;
//            }
//            if(paras.size()==2){
//                T first=paras[0]==""?0
//            }
//            if(paras.size()==3){
//
//            }
//        }
	}
	/**\brief encapsule some IO operation
	  */
	namespace io {
        /**\brief check file for existence
          */
        bool fileExists(std::string fileName);
        /**\brief print console progress bar(better print a endl first before progress bar started)
          *\param[in] value the current progress value(0~100)
          *\param[in] delimiterL the left delimiter
          *\param[in] increment the progress symbol
          *\param[in] delimiterR the right deimiter
          *\param[in] status a string to declare current status(eg. "processing","reading",etc)
          *\note: when value is 100, it will automatically go to next line;
          *\progress update within 0-100 will update in the same line
          *\note the update frequency should not be very
          *\for(int i=0;i<total;++i){
          *\	if(i%(total/100)==0)lyxutils::consoleProgressBar(i*100/total);
          *\	//some enduring repetitive procedure
          *\}
          */
        void consoleProgressBar(int value, std::string delimiterL = "[", std::string increment = ">",
                                std::string delimiterR = "]", std::string status = "");
        /**
         * print console progress bar,this version is suggested to replace the older version:
         * void consoleProgressBar(int,std::string,std::string,std::string,std::string)
         * because this one can set progress bar width and automatically avoid redundant update
         * @param value current state,value ranges from 1-total
         * @param total total number of loops
         * @param bar_width progress bar width(default is 50),left/right delimiter and xx% not included
         * @param delimiterL left delimiter(default is '[')
         * @param increment increment char(default is '>')
         * @param space white space char(default is ' ')
         * @param delimiterR right delimiter(default is ']')
         */
        void consoleProgressBar(int value, int total, int bar_width=50, char delimiterL='[',
                                char increment='>', char space=' ', char delimiterR=']');
        /**\brief an easy wrapper function used to print a vector(the vector element type need to implement << operator)
         * @tparam T type of vector element
         * @tparam alloc allocator type
         * @param vec output vector
         * @param sep separator
         * @param tight use tight layout(do not add space to force alignnment)
         * @param w width of each place
         * @param width total display column
         */
        template<typename T,typename alloc=std::allocator<T> >
        void printVector(const std::vector<T,alloc> &vec,std::string sep=",", bool tight=true, int w=6, int width=50) {
            if(vec.size()==0){
                std::cout<<"[]"<<std::endl;
            }
            std::cout.setf(std::ios::right);
            if(vec.size()<=width) {
                for_each(vec.begin(), vec.end(), [&](const T &t) { if(!tight)std::cout.width(w);std::cout << t << sep; });
            }
            else{
                int left=width/2,right=width-left;
                for(int i=0;i<left;++i){
                    if(!tight)std::cout.width(w);
                    std::cout<<vec[i]<<sep;
                }
                std::cout<<str_utils::center(".",w,'.')<<sep;
                for(int i=vec.size()-right;i<vec.size();++i){
                    if(!tight)std::cout.width(w);
                    std::cout<<vec[i]<<sep;
                }
            }
            std::cout<<std::endl;
        }
        template<typename T,typename alloc=std::allocator<T> >
        void print2DVector(const std::vector<std::vector<T,alloc> > &vec, std::string sep=",", bool tight=true, int w=6, int width=14, int height=30) {
            if(vec.size()==0 || vec[0].size()==0){
                std::cout<<"[[]]"<<std::endl;
            }
            if(vec.size()<=height){
                std::for_each(vec.begin(),vec.end(),[&](const std::vector<T,alloc> &v){printVector(v,sep,tight,w,width);});
            }
            else{
                height=height-3;//3 lines for dots
                int top=height/2,bottom=height-top;
                for(int i=0;i<top;++i){
                    printVector(vec[i],sep,tight,w,width);
                }
                int realWidth=width>vec[0].size()?vec[0].size():width;
                int total=w*realWidth+realWidth-1;
                int left=total/3,right=total/3,middle=total-left-right;
                if(!tight){
                    std::string line=str_utils::left(str_utils::multiply(" ",w/2)+".",left,' ')+str_utils::center(".",middle,' ')+str_utils::right("."+str_utils::multiply(" ",w/2),right,' ');
                    for(int i=0;i<3;++i)std::cout<<line<<std::endl;
                }
                else{
                    for(int i=0;i<3;++i)std::cout<<"\t.\t.\t."<<std::endl;
                }
                for(int i=vec.size()-bottom;i<vec.size();++i){
                    printVector(vec[i],sep,tight,w,width);
                }
            }
        }
        /**\brief create a folder named folderPath if it does not exist
          *\param[in] folderPath the folder path to be created
          *\return true if creation succeeded or the folder already exists;false otherwise
          */
        bool createFolder(const std::string folderPath);

        /**\brief open a log file for writing log
          *\usage: std::ofstream *log = dbg::openLogFile(file);
          *\       dbg::writeLog(*log,message);
          */
        std::ofstream *openLogFile(std::string fileName);

        /**\brief write message to opened log file
          *\param[in] log the file stream to write to
          *\param[in] message the message to write
          *\param[in] end the message end(default end is "\n",automatically carriage return)
          */
        void writeLog(std::ofstream &log, std::string message, std::string end = "\n");

        /**\brief read csv file
		  *\param[in] fileName the input file name
		  *\param[out] table the 2d vector to store the data.table is cleared before reading
		  *\param[in] sep the separator to separate each column, usually ","
		  *\param[out] report the reading process report(record open file name,time,bad lines, etc)
		  *\param[in] fields the number of fields to read(-1 means use total number of fields in the first line of the file)
		  *\param[in] read_header whether to read first line as header or not
		  *\param[out] headers the headers of the table(if read_header is true)
		  *\param[in] numberOfRedirect number of fields to be redirected(not stored in table),-1 for all fields
		  *\param[in] directFunc the function to redirect first numberOfRedirect of fields
		  *\note:parameter of directFunc must have the same dimension as numberOfRedirect
		  *\the last parameter directFuncs can be defined as follows if trying to redirect to point cloud(cloud_in is a global pc pointer)
			void direct(const std::vector<float> &vec) {
				pcl::PointXYZRGB pt = pcl::PointXYZRGB();
				pt.x = vec[0];
				pt.y = vec[1];
		    	pt.z = vec[2];
				pt.r = (uint8_t)(vec[3]);
				pt.g = (uint8_t)(vec[4]);
				pt.b = (uint8_t)(vec[5]);
				cloud_in->points.push_back(pt);
			}
		*\addition: if the columns of data is less than required, it will fill with 0s;
		*\			if data column is larger than required, the rest will be abbreviated
        *\          this function does not deal with missing value, e.g.if 4 data a line
         *          dataline:2.3,,5.7,8.9=>can not discover missing value at second place, will read as 2.3,5.7,8.9,0
		*/
        bool __read_csv(const std::string &fileName, std::vector<std::vector<float> > &table,
                const std::string &sep, bool report, std::ostream &ssReport, int fields,
                bool read_header, std::vector<std::string> &headers, int numberOfRedirect,
                std::function<void(const std::vector<float> &row)> directFunc);
        /**\brief read csv file to 2d vector,separated by sep
         */
        bool read_csv(std::string fileName, std::vector<std::vector<float> > &table, const std::string &sep = ",");
        bool read_csv(std::string fileName, std::vector<std::vector<float> > &table, const std::string &sep, std::ostream &out);
        /**\brief read csv file to 2d vector,separated by sep. take first line as header
         */
        bool read_csv(std::string fileName, std::vector<std::vector<float> > &table, std::vector<std::string> &header,const std::string &sep = ",");
        bool read_csv(std::string fileName, std::vector<std::vector<float> > &table, std::vector<std::string> &header,const std::string &sep,std::ostream &out);
        /**\brief read csv file if only data line satisfies a certain condition
         * note:pred is a lambda with signature [](const std::vector<float> &row)->bool{...} (recommended)
         */
        template<class Pred>
        bool read_csv_if(std::string fileName,std::vector<std::vector<float> > &table,Pred pred,const std::string sep=","){
            auto read_line=[&](const std::vector<float> &row)->void{
                if(row.size()>0&&pred(row)){
                    table.push_back(std::vector<float>());
                    for(auto it=row.begin();it!=row.end();++it){
                        (table.end()-1)->push_back(*it);
                    }
                }
            };
            std::vector<std::vector<float> > temp;
            std::vector<std::string> header;
            return __read_csv(fileName,temp,sep,false,std::cout,-1,false,header,-1,read_line);
        }
        template<class Pred>
        bool read_csv_if(std::string fileName,std::vector<std::vector<float> > &table,Pred pred,const std::string sep, std::ostream &out){
            auto read_line=[&](const std::vector<float> &row)->void{
                if(row.size()>0&&pred(row)){
                    table.push_back(std::vector<float>());
                    for(auto it=row.begin();it!=row.end();++it){
                        (table.end()-1)->push_back(*it);
                    }
                }
            };
            std::vector<std::vector<float> > temp;
            std::vector<std::string> header;
            return __read_csv(fileName,temp,sep,true,out,-1,false,header,-1,read_line);
        }
        /**
         * read csv file and process each row using given function,this enables customized reading process
         * @param fileName csv file name(without header row)
         * @param procRow process function-can be a normal function,functor and lambda expression
         * @param sep separator sign
         * @return
         */
        bool read_csv_customize(std::string fileName,std::function<void(const std::vector<float>&)> procRow,std::string sep=",");
        /**
         * read csv file and process each row using given function,this enables customized reading process
         * @param fileName csv file name(without header row)
         * @param procRow process function-can be a normal function,functor and lambda expression
         * @param sep separator sign
         * @param out log output when reading csv file(can be a ofstream or std::cout,etc)
         * @return
         */
        bool read_csv_customize(std::string fileName,std::function<void(const std::vector<float>&)> procRow,std::string sep,std::ostream &out);
        /**\brief read point cloud from csv file(without header)
         */
        template<class PointXYZ, class allocator=std::allocator<PointXYZ> >
        bool read_pts(std::string fileName, std::vector<PointXYZ, allocator> &pts, const std::string &sep = ",") {
            auto toPts = [&pts](const std::vector<float> &row) -> void {
                PointXYZ pt;
                pt.x = row[0];
                pt.y = row[1];
                pt.z = row[2];
                pts.push_back(pt);
            };
            std::vector<std::vector<float> > table;
            std::vector<std::string> header;
            pts.clear();
            return __read_csv(fileName, table, sep, false, std::cout, -1, false, header, 3, toPts);
        }
        template<class PointXYZ, class allocator=std::allocator<PointXYZ> >
        bool read_pts(std::string fileName, std::vector<PointXYZ, allocator> &pts, const std::string &sep, std::ostream &out){
            auto toPts = [&pts](const std::vector<float> &row) -> void {
                PointXYZ pt;
                pt.x = row[0];
                pt.y = row[1];
                pt.z = row[2];
                pts.push_back(pt);
            };
            std::vector<std::vector<float> > table;
            std::vector<std::string> header;
            pts.clear();
            return __read_csv(fileName, table, sep, true, out, -1, false, header, 3, toPts);
        }
        /**\brief read point cloud and fields from csv file(without header)
         */
        template<class PointXYZ, class allocator=std::allocator<PointXYZ> >
        bool read_pts_field(std::string fileName, std::vector<PointXYZ, allocator> &pts,
                                   std::vector<std::vector<float> > &fields, const std::string &sep = ",") {
            auto toPts = [&pts](const std::vector<float> &row) -> void {
                PointXYZ pt;
                pt.x = row[0];
                pt.y = row[1];
                pt.z = row[2];
                pts.push_back(pt);
            };
            std::vector<std::string> header;
            pts.clear();
            return __read_csv(fileName, fields, sep, false, std::cout, -1, false, header, 3, toPts);
        }
        template<class PointXYZ, class allocator=std::allocator<PointXYZ> >
        bool read_pts_field(std::string fileName, std::vector<PointXYZ, allocator> &pts,
                std::vector<std::vector<float> > &fields, const std::string &sep,std::ostream &out){
            auto toPts = [&pts](const std::vector<float> &row) -> void {
                PointXYZ pt;
                pt.x = row[0];
                pt.y = row[1];
                pt.z = row[2];
                pts.push_back(pt);
            };
            std::vector<std::string> header;
            pts.clear();
            return __read_csv(fileName, fields, sep, true, out, -1, false, header, 3, toPts);
        }
        /**\brief read point cloud and fields from csv file(with header)
         */
        template<class PointXYZ, class allocator=std::allocator<PointXYZ> >
        bool read_pts_field(std::string fileName, std::vector<PointXYZ, allocator> &pts, std::vector<std::string> &header,
                       std::vector<std::vector<float> > &fields, const std::string &sep = ",") {
            auto toPts = [&pts](const std::vector<float> &row) -> void {
                PointXYZ pt;
                pt.x = row[0];
                pt.y = row[1];
                pt.z = row[2];
                pts.push_back(pt);
            };
            pts.clear();
            return __read_csv(fileName, fields, sep, false,std::cout, -1, true, header, 3, toPts);
        }
        template<class PointXYZ, class allocator=std::allocator<PointXYZ> >
        bool read_pts_field(std::string fileName, std::vector<PointXYZ, allocator> &pts, std::vector<std::string> &header,
                            std::vector<std::vector<float> > &fields, const std::string &sep,std::ostream &out){
            auto toPts = [&pts](const std::vector<float> &row) -> void {
                PointXYZ pt;
                pt.x = row[0];
                pt.y = row[1];
                pt.z = row[2];
                pts.push_back(pt);
            };
            pts.clear();
            return __read_csv(fileName, fields, sep, true, out, -1, true, header, 3, toPts);
        }
        /**\brief read colored point cloud from csv file(without header)
         */
        template<class PointXYZRGB, class allocator=std::allocator<PointXYZRGB> >
        bool read_ptsrgb(std::string fileName, std::vector<PointXYZRGB, allocator> &pts, const std::string &sep = ",") {
            auto toPts = [&pts](const std::vector<float> &row) -> void {
                PointXYZRGB pt;
                pt.x = row[0];
                pt.y = row[1];
                pt.z = row[2];
                pt.r = row[3];
                pt.g = row[4];
                pt.b = row[5];
                pts.push_back(pt);
            };
            std::vector<std::vector<float> > table;
            std::vector<std::string> header;
            pts.clear();
            return __read_csv(fileName, table, sep, false, std::cout, -1, false, header, 6, toPts);
        }
        template<class PointXYZRGB, class allocator=std::allocator<PointXYZRGB> >
        bool read_ptsrgb(std::string fileName, std::vector<PointXYZRGB, allocator> &pts, const std::string &sep,std::ostream &out){
            auto toPts = [&pts](const std::vector<float> &row) -> void {
                PointXYZRGB pt;
                pt.x = row[0];
                pt.y = row[1];
                pt.z = row[2];
                pt.r = row[3];
                pt.g = row[4];
                pt.b = row[5];
                pts.push_back(pt);
            };
            std::vector<std::vector<float> > table;
            std::vector<std::string> header;
            pts.clear();
            return __read_csv(fileName, table, sep, true, out, -1, false, header, 6, toPts);
        }
        /**\brief read colored point cloud and fields from csv file(without header)
         */
        template<class PointXYZRGB, class allocator=std::allocator<PointXYZRGB> >
        bool read_ptsrgb_field(std::string fileName, std::vector<PointXYZRGB, allocator> &pts,
                                      std::vector<std::vector<float> > &fields, const std::string &sep = ",") {
            auto toPts = [&pts](const std::vector<float> &row) -> void {
                PointXYZRGB pt;
                pt.x = row[0];
                pt.y = row[1];
                pt.z = row[2];
                pt.r = row[3];
                pt.g = row[4];
                pt.b = row[5];
                pts.push_back(pt);
            };
            std::vector<std::string> header;
            pts.clear();
            return __read_csv(fileName, fields, sep, false, std::cout, -1, false, header, 6, toPts);
        }
        template<class PointXYZRGB, class allocator=std::allocator<PointXYZRGB> >
        bool read_ptsrgb_field(std::string fileName, std::vector<PointXYZRGB, allocator> &pts,
                               std::vector<std::vector<float> > &fields, const std::string &sep,std::ostream &out){
            auto toPts = [&pts](const std::vector<float> &row) -> void {
                PointXYZRGB pt;
                pt.x = row[0];
                pt.y = row[1];
                pt.z = row[2];
                pt.r = row[3];
                pt.g = row[4];
                pt.b = row[5];
                pts.push_back(pt);
            };
            std::vector<std::string> header;
            pts.clear();
            return __read_csv(fileName, fields, sep, true, out, -1, false, header, 6, toPts);
        }
        /**\brief read colored point cloud and fields from csv file(with header)
         */
        template<class PointXYZRGB, class allocator=std::allocator<PointXYZRGB> >
        bool read_ptsrgb_field(std::string fileName, std::vector<PointXYZRGB, allocator> &pts,
                                      std::vector<std::string> &header,
                                      std::vector<std::vector<float> > &fields, const std::string &sep = ","){
            auto toPts = [&pts](const std::vector<float> &row) -> void {
                PointXYZRGB pt;
                pt.x = row[0];
                pt.y = row[1];
                pt.z = row[2];
                pt.r = row[3];
                pt.g = row[4];
                pt.b = row[5];
                pts.push_back(pt);
            };
            pts.clear();
            return __read_csv(fileName, fields, sep, false, std::cout, -1, true, header, 6, toPts);
        }
        template<class PointXYZRGB, class allocator=std::allocator<PointXYZRGB> >
        bool read_ptsrgb_field(std::string fileName, std::vector<PointXYZRGB, allocator> &pts,
                               std::vector<std::string> &header,
                               std::vector<std::vector<float> > &fields, const std::string &sep,std::ostream &out){
            auto toPts = [&pts](const std::vector<float> &row) -> void {
                PointXYZRGB pt;
                pt.x = row[0];
                pt.y = row[1];
                pt.z = row[2];
                pt.r = row[3];
                pt.g = row[4];
                pt.b = row[5];
                pts.push_back(pt);
            };
            pts.clear();
            return __read_csv(fileName, fields, sep, true, out, -1, true, header, 6, toPts);
        }
        /**\brief write csv file
          */
        template<class DataType>
        void write_csv(std::string fileName, const std::vector<std::vector<DataType>> &table, std::string sep, bool show_progress=false) {
            std::ofstream ofs(fileName);
            ofs << std::setprecision(8);
            for(int i=0;i<table.size();++i){
                ofs<<(*table[i].begin());
                std::for_each(table[i].begin()+1,table[i].end(),[&](const DataType &data){ofs<<sep<<data;});
                ofs<<std::endl;
                if(show_progress)consoleProgressBar(i+1,table.size());
            }
            if (ofs.is_open()) {
                ofs.close();
            }
        }
        /**\brief write points to csv file
          */
        template<class PointXY, class allocator=std::allocator<PointXY> >
        void points_to_csv(std::string fileName, const std::vector<PointXY, allocator> &points,
                std::string sep = ",", int precision = 8, bool show_progress=false) {
            std::ofstream ofs(fileName);
            ofs << std::setprecision(precision);
            for (int i = 0; i < points.size(); ++i) {
                ofs << points[i].x << sep << points[i].y << std::endl;
                if(show_progress)consoleProgressBar(i+1,points.size());
            }
            if (ofs.is_open()) {
                ofs.close();
            }
        }
        /**\brief write 3D points to file
          */
        template<class PointXYZ, class allocator = std::allocator<PointXYZ> >
        void points3D_to_csv(std::string fileName, const std::vector<PointXYZ, allocator> &points,
                std::string sep = ",", int precision = 8, bool show_progress=false) {
            std::ofstream ofs(fileName);
            ofs << std::setprecision(precision);
            for (int i = 0; i < points.size(); ++i) {
                ofs << points[i].x << sep << points[i].y << sep << points[i].z << std::endl;
                if(show_progress)consoleProgressBar(i+1,points.size());
            }
            if (ofs.is_open()) {
                ofs.close();
            }
        }
        /**\brief write 3D colored points to file
        */
        template<class PointXYZRGB, class allocator = std::allocator<PointXYZRGB> >
        void points3DRGB_to_csv(std::string fileName, std::vector<PointXYZRGB, allocator> &points,
                std::string sep = ",", int precision = 8, bool show_progress=false) {
            std::ofstream ofs(fileName);
            ofs << std::setprecision(precision);
            for (int i = 0; i < points.size(); ++i) {
                ofs << points[i].x << sep << points[i].y << sep << points[i].z << sep
                    << points[i].r << sep << points[i].g << sep << points[i].b << std::endl;
                if(show_progress)consoleProgressBar(i+1,points.size());
            }
            if (ofs.is_open()) {
                ofs.close();
            }
        }
        /**\brief command line option
         * define long option name, whether it has argument and what is the corresponding short option
         * optName-define the long option name,i.e.help
         * hasArg -define whether the long option requires argument:0-no argument,1-optional argument,2-require argument
         * abbr   -define the corresponding short version of the option:NULL-no short version,not NULL-point to the short option
         * i.e: CLOption opt={"target-location",2,new char('t')};
         */
        struct CLOption {
            std::string optName;
            int hasArg;
            bool hasAbbr;
            char abbr;
        };
        /**\brief command line parser
         * \usage:
         *      //1.if no option pattern is given,every option requires an argument or optional argument
         *      CLParser parser;
         *      std::string cmd="cmd --target a.txt -v --when=always b.txt";
         *      try{
         *          parser.parse(cmd);
         *          parser.hasOption("v");//=>bool:true
         *          parser.getOptionArg("when");//=>string:"always"
         *          parser.getParameters();//=>vector of strings,contains "b.txt",
         *          parser.getAllOptions();//=>map of option-argument pair
         *      }catch(std::invalid_argument ia){cout<<ia.what()<<endl;}
         *      //2.if option pattern is given,all applicable short option will be replaced by longer ones
         *      CLOption options={
         *          {"verbose",0,new char('v')},//--verbose is equivalent to -v, this option requires no argument
         *          {"when",1,NULL},            //--when has no abbreviated form, --when=arg to assign argument to option
         *          {"target",2,new char('t')}  //--target is equivalent to -t, this option requires an argument
         *      };
         *      try{...}catch(std::invalid_argument ia){cout<<ia.what()<<endl;}
         *      //3.for sub commands parse:
         *      if(argc>1&&string(argv[1])=="sub1"){dealSub1(argc-1,argv[1:]);}
         *      else if(argc>1&&string(argv[1])=="sub2"){dealSub2(argc-1,argv[1:]);}
         *      ...
         *      else dealCommand(argc,argv);
         */
        class CLParser {
        public:
            CLParser();

            /**\brief parse a line of command
             * split command line using space or tab and feed to function parse(const vector<string> &argv) to parse
             * this function throws exceptions when:
             *      1.unrecognized option(only when option pattern is given)->invalid_argument
             *      2.option has no argument(only when argument is required)->invalid_argument
             */
            void parse(std::string command);

            /**\brief parse command line
             * this function throws exceptions when:
             *      1.unrecognized option(only when option pattern is given)->invalid_argument
             *      2.option has no argument(only when argument is required)->invalid_argument
             */
            void parse(int argc, char **argv);

            /**\brief concatenate string of `commands` and string from file `filename` and parse
			 * \string from "filename" will be preprocessed to force ":" equivalent to "=",allows comments(by "//" or "#") 
			 * and omit line not starting with "--";in parsing command from file mode,only option-arg pairs are recommended
			 * 
             * this function throws exceptions when:
             *      1.unrecognized option(only when option pattern is given)->invalid_argument
             *      2.option has no argument(only when argument is required)->invalid_argument
             *      3.file not exist->runtime_error
             *      4.open file failed->runtime_error
             */
            void parse(std::string commands, std::string filename);

            /**\brief set options and arguments pattern
             */
            void setOptionPattern(int num_opts, const lyxutils::io::CLOption *options);

            void clearOptionPattern();

            /**\brief check existence of option in the given command
             */
            bool hasOption(std::string optName);

            /**\brief get the option value
             * this option does not throw exception,user must use hasOption first to check the existence of optName
             * if given option doesn't exist,return ""
             */
            std::string getOptionArg(std::string optName);

            /**\brief get all options parsed
             */
            std::map<std::string, std::string> getAllOptions();

            /**\brief get all parameters
             */
            std::vector<std::string> getParameters();

            /**\brief add example of usages to help message
             */
            void addUsageExample(const std::string &example);
            /**\brief generate command line help message
             * @param basicUsage the current command of this parser
             * @param requireObj whether this command required an object
             * @param description the description of this command
             * @param _optNames all the available option names(if option pattern is set,the number of names should equal
             *                 to the number of predefined options, or a logic_error will be throw; Each of the option
             *                 name should also equal the predefined option name, or a logic_error will be throw
             * @param optDescription description for each option
             * @param optArgName argument name of each option
             * @return help message string
             */
            std::string generateHelp(std::string basicUsage, std::string description,
                                     const std::vector<std::string> &_optNames,
                                     const std::vector<std::string> &optDescription,
                                     const std::vector<std::string> &optArgName);
            /**\brief generate command line help message(using help message set by last time)
             * @return help message string
             */
            std::string generateHelp();

        private:
            bool _pattern_set = false;
            std::map<std::string, std::string> _opt_val;
            std::vector<std::string> _parameters;
            std::map<std::string, std::string> _short2verbose;
            std::map<std::string, CLOption> _opt_patterns;
            std::vector<std::string> _usageExamples;
            //last time help message record
            std::string _last_usage;
            std::string _last_description;
            std::vector<std::string> _last_optNames;
            std::vector<std::string> _last_optDescription;
            std::vector<std::string> _last_argNames;

            void parse(std::vector<std::string> argv);
            void split_command(std::string command,std::vector<std::string> &argv);
        };
    }
#ifdef enable_deprecated
	/**\brief encapsulate some alorithms and Data structures
	  */
	namespace algorithms 
	{
		const double PI = 3.1415926535897932;

		/**\brief compute the average x y z coordinates of points contained in pList
		  *\param[in] pList the vector contains the points
		  *\param[out] average_x,average_y,average_z the computed averages
		  */
		template<class PointXYZ,class allocator=std::allocator<PointXYZ> >
		void computeAverages(const std::vector<PointXYZ,allocator> &pList, double &average_x, double &average_y, double &average_z)
		{
			average_x = 0;
			average_y = 0;
			average_z = 0;
			for (int i = 0; i < pList.size(); ++i) {
				average_x += pList[i].x;
				average_y += pList[i].y;
				average_z += pList[i].z;
			}
			if (pList.size() > 0) {
				average_x = average_x / pList.size();
				average_y = average_y / pList.size();
				average_z = average_z / pList.size();
			}
		}
		/**\brief find the distance of point p to plane defined by ax+by+cz+d=0, where a^2+b^2+c^2=1
		  */
		template<class PointXYZ>
		double distanceToPlane(const PointXYZ &p, double a, double b, double c, double d)
		{
			return abs(a*p.x + b*p.y + c*p.z + d);
		}
		/**\brief find the distance of point (x,y,z) to plane defined by ax+by+cz+d=0
		  *\where a^2+b^2+c^2=1
		  */
		template<typename T>
		T distanceToPlane(T x,T y,T z,T a,T b,T c,T d){
            return abs(a*x + b*y + c*z + d);
		}
		/**\brief find the average distance of given points to a normalized plane with given parameters
		  *\param[in] pList the point vector
		  *\param[in] a,b,c,d the plane coefficient
		  *\param[out] average_dis the average distance of points to plane
		  */
		template<class PointXYZ, class allocator = std::allocator<PointXYZ> >
		void findAverageDistanceToPlane(const std::vector<PointXYZ, allocator> &pList, const double &a,
			const double &b, const double &c, const double &d, double &average_dis)
		{
			average_dis = 0;
			for (int i = 0; i < pList.size(); ++i) {
				average_dis += distanceToPlane(pList[i].x, pList[i].y, pList[i].z, a, b, c, d);
			}
			average_dis = average_dis / pList.size();
		}
		/**\brief region growing algorithm,give a seed unit and grow out a whole region.mark each unit being growned
		  *\param[in] seed the given seed point(type is unit,define whichever type you think can represent the growing unit)
		  *\param[in] findNeighbors a function or functor that can find the neighbors of a given unit
		  *\param[in] growingCondition the growing condition for this algorithm
		  *\param[out] growingResult the growing result,a vector that contains all the growing resultant units
		  *\function or functor declarations:
		  *\std::vector<Unit>* findNeighbors(const Unit &center)
		  *\bool growingCondition(const Unit &center,const Unit &neighbor)
		  *\
		  *\note: this tenmplate function is ill-designed, UnaryOperation and BinaryOperation should respectively 
		  *\		be functions as:`std::vector<Unit>* ()(Unit &center)` and `bool ()(Unit &center,Unit &target)`
		  */
		template<class Unit, class UnaryOperation, class BinaryPredict>
		void regionGrowing(Unit &seed, UnaryOperation findNeighbors,
			BinaryPredict growingCondition, std::vector<Unit> &growingResult)
		{
			std::vector<Unit> *recursion = new std::vector<Unit>();
			growingResult.push_back(seed);
			recursion->push_back(seed);
			while (recursion->size() > 0)
			{
				Unit indexesTemp = (*recursion)[recursion->size() - 1];
				recursion->pop_back();
				std::vector<Unit>* neighbors = findNeighbors(indexesTemp);
				for (auto iterator = neighbors->begin(); iterator != neighbors->end(); iterator++) {
					if (growingCondition(indexesTemp, *iterator)) {
						growingResult.push_back(*iterator);
						recursion->push_back(*iterator);
					}
				}
				delete neighbors;
			}
			delete recursion;
		}
        /**\brief find the orientation of the ordered triplet (p1,p2,p3)
          *\function returns following values:
          *\0-->if p1,p2,p3 are colinear (or if at least two of them are indentical)
          *\1-->if clockwise
          *\2-->if counterclockwise
          */
        template<class PointXY>
        int orientation(const PointXY &p1, const PointXY &p2, const PointXY &p3)
        {
            double val = (p2.y - p1.y) * (p3.x - p2.x) - (p2.x - p1.x) * (p3.y - p2.y);
            if (val == 0) return 0;  // colinear
            return (val > 0) ? 1 : 2; // clock or counterclock wise
        }
        /**\brief compare the two polar angles x-p1-p2 and x-p1-p3(counterclockwise)
          *\return true if x-p1-p2 is smaller than x-p1-p3(or if p1p2<p1p3 when p1,p2,p3 are colinear)
          */
        template<class PointXY>
        bool compareAngle(const PointXY &p1, const PointXY &p2, const PointXY &p3)
        {
            int o = orientation(p1, p2, p3);
            if (o == 0) {
                //return distance2DSquared(p1, p2) < distance2DSquared(p1, p3);
                //return distance2D(p1, p2) < distance2D(p1, p3);
                //assume p1 is the left most one
                double v12_x = p2.x - p1.x;
                double v12_y = p2.y - p1.y;
                double v23_x = p3.x - p2.x;
                double v23_y = p3.y - p2.y;
                double innerProduct = v12_x*v23_x + v12_y*v23_y;
                return innerProduct > 0;
            }
            return (o == 2) ? true : false;
        }
		/**\brief use Graham Scan to get the convex hull(anti-clockwise starting from min x) border of a given set of points
		  *\param[in] pList the given set of points
		  *\param[out] border the output set of border points
		  */
		template<class PointXY, class allocator=std::allocator<PointXY> >
		void findConvexHullBorder(const std::vector<PointXY,allocator> &pList, std::vector<PointXY,allocator> &border)
		{
			if (pList.size() < 3)
				throw std::runtime_error("number of points too small(<3) to define a convex hull");
			//copy to border,operate on border
			double x_min = pList[0].x;
			int x_min_index = 0;
			border.push_back(pList[0]);
			for (int i = 1; i < pList.size(); ++i) {
				border.push_back(pList[i]);
				if (pList[i].x < x_min || pList[i].x == x_min&&pList[i].y < pList[x_min_index].y) {
					x_min = pList[i].x;
					x_min_index = i;
				}
			}
			//move smallest x to head
			PointXY pointTemp = border[x_min_index];
			border[x_min_index] = border[0];
			//order anti-clock-wise
			//sort(border.begin() + 1, border.end(), [&pointTemp](PointType *p1, PointType *p2)
			//->bool {return compareAngle(*pointTemp, *p1, *p2); });
			border.erase(border.begin());
			merge_sort(border, [&pointTemp](PointXY p1, PointXY p2)
				->bool {return compareAngle(pointTemp, p1, p2); });
			border.insert(border.begin(), pointTemp);
			//remove collinear point, maintain only the furthest
			for (auto it = border.begin() + 1; it != border.end();) {
				if ((it + 1) != border.end() && orientation(pointTemp, *it, *(it + 1)) == 0) {
					it = border.erase(it);
				}
				else {
					it++;
				}
			}
			if (border.size() < 3)
				throw std::runtime_error("number of border points too small(<3)");
			if (border.size() == 3)
			{
				return;
			}
			//remove concave vertex
			auto first = border.begin() + 1;
			auto second = first + 1;
			auto third = second + 1;
			while (third != border.end()) {
				int orien = orientation(*first, *second, *third);
				if (orien != 2) {
					border.erase(second);
					second = first;
					first--;
					third = second + 1;
				}
				else {
					first++;
					second++;
					third++;
				}
			}
		}
		/**\brief find convex hull using javis march algorithm, the border point is in anti-clockwise order starting from mkinimum y value point
		  *\param[in] pList input point vector
		  *\param[out] border output vertices
		  *\note: the point class should have default constructor that requires no argument,should have clone & copy constructor
		  *\	  the temporal complexity is O(nH), where H is the number of vertices;
		  *\	  should be faster than graham-O(nlogn) when number of vertices are limited
		  */
		template<class PointXY, class allocator = std::allocator<PointXY> >
		void javisMarch(const std::vector<PointXY, allocator> &pList, std::vector<PointXY, allocator> &border)
		{
			if (pList.size() < 3)
				throw std::runtime_error("number of points too small(<3) to define a convex hull");
			//get the vector comprised by two points
			auto vector2P = [](const PointXY &p1, const PointXY &p2)->PointXY {
				PointXY result;
				result.x = p2.x - p1.x;
				result.y = p2.y - p1.y;
				return result;
			};
			auto turningAngle = [&vector2P](const PointXY *p1, const PointXY *p2, const PointXY *p3)->float
			{
				PointXY v1;
				PointXY v2;
				if (p1 == NULL) {
					v1.x = 1;
					v1.y = 0;
				}
				else {
					v1 = vector2P(*p1, *p2);
				}
				v2 = vector2P(*p2, *p3);
				return (v1.x*v2.x + v1.y*v2.y) / (sqrt(v1.x*v1.x + v1.y*v1.y)*sqrt(v2.x*v2.x + v2.y*v2.y));
			};
			double yMin = pList[0].y;
			int min_index = 0;
			std::vector<const PointXY*> remainPts;
			for (int i = 0; i < pList.size(); ++i) {
				remainPts.push_back(&pList[i]);
				if (pList[i].y < yMin) {
					yMin = pList[i].y;
					min_index = i;
				}
			}
			const PointXY *pTemp = remainPts[0];
			remainPts[0] = remainPts[min_index];
			remainPts[min_index] = pTemp;
			const PointXY *p1, *p2;
			p1 = NULL;
			p2 = remainPts[0];
			border.push_back(*p2);
			const PointXY *starter = p2;
			while (true) {
				if (remainPts.size() < 2)break;
				float maxCos = -2.0;
				auto p3_pntr = remainPts.begin();
				for (auto it = remainPts.begin(); it != remainPts.end(); ) {
					if ((*it)->x == p2->x && (*it)->y == p2->y) {
						if (it != remainPts.begin()) {
							it = remainPts.erase(it);
						}
						else {
							++it;
						}
					}
					else {
						float cosVal = turningAngle(p1, p2, *it);
						if (cosVal > maxCos) {
							maxCos = cosVal;
							p3_pntr = it;
						}
						++it;
					}
				}
				p1 = p2;
				p2 = *p3_pntr;
				border.push_back(*p2);
				remainPts.erase(p3_pntr);
				if (p2==starter) {
					border.pop_back();
					break;
				}
			}
		}
		/**\brief creationg mask image that marks the interior points inside a convex hull border
		  *\param[in] border the convex hull border (the point coordinate x y should be integer
		  *\param[in] rows,cols the rows and cols of the mask image
		  *\param[out] access the function that sets image value at x,y position(like access(x,y)=data),
		  *\	eg. auto accessMat=[&mask](int x,int y)->uchar&{return mask.at<uchar>(x,y);};
		  *\		//above lambda function can used to access a cv::Mat of type CV_8UC1
		  *\		convexhullIntereriorMask(border,mask.rows,mask.cols,accessMat,255);
		  *\param[in] value the foreground value(the value of points inside convex hull)
		  */
		template<class PointXY,class allocator=std::allocator<PointXY>,class DataType,class Accessor>
		void convexhullIntereriorMask(const std::vector<PointXY, allocator> &border, int rows, int cols, Accessor access, DataType value)
		{
			std::vector<int> left(rows, cols), right(rows, -1);//record left border y value and right border y value
			int xMinIdx = 0;
			int xMin = border[0].x;
			for (int i = 1; i < border.size();++i) {
				if (border[i].x < xMin) {
					xMin = border[i].x;
					xMinIdx = i;
				}
			}
			for (int i = xMinIdx; i < xMinIdx + border.size(); ++i) {//start from min x to max x then back to min x
				int x1 = border[i%border.size()].x;
				int y1 = border[i%border.size()].y;
				int x2 = border[(i + 1) % border.size()].x;
				int y2 = border[(i + 1) % border.size()].y;
				if (x1 < x2) {//left side
					for (int k = x1; k < x2; ++k) {
						left[k] = float(y1*(x2 - k) + y2*(k - x1)) / (x2 - x1);
					}
				}
				else if (x1 == x2) {
					left[x1] = y1;
					right[x2] = y2;
				}
				else {
					for (int k = x1; k > x2; --k) {
						right[k] = float(y1*(k - x2) + y2*(x1 - k)) / (x1 - x2);
					}
				}
			}
			for (int i = 0; i < rows; ++i) {//fill the mask image
				for (int j = left[i]; j < right[i]; ++j) {
					access(i, j) = value;
				}
			}
		}
		/**\brief calculate the area of a convex hull
		  *\param[in] border the border p[oint list��can be derived from lyxutils::algorithm::findConvexHull
		  *\note:the point list is in the order of clockwise or anti-clockwise
		  */
		template<class PointXY,class allocator=std::allocator<PointXY> >
		double getConvexHullArea(const std::vector<PointXY,allocator> &border)
		{
			if (border.size() < 3)
				throw std::runtime_error("number of points too small(<3) to enclose");
			auto assign = [](const PointXY &from, PointXY &to) {to.x = from.x; to.y = from.y; };//�����û�п�¡���캯��
			PointXY p1, p2, p3;
			double area = 0;
			int i = 1;
			assign(border[0], p1);
			while (i + 1 < border.size()) {
				assign(border[i], p2);
				assign(border[i + 1], p3);
				area += triangleArea(p1, p2, p3);
				++i;
			}
			return area;
		}
		/**\brief calculate area of a triangle(2D)
		*/
		template<class PointXY>
		double triangleArea(const PointXY &vertice1, const PointXY &vertice2, const PointXY &vertice3)
		{
			double m_x = vertice2.x - vertice1.x;
			double m_y = vertice2.y - vertice1.y;
			double n_x = vertice3.x - vertice1.x;
			double n_y = vertice3.y - vertice1.y;
			return abs(m_x*n_y - m_y*n_x) / 2;
		}
		/**\brief find the squared Euclidean distance between two points
		*/
		template<class PointXY>
		double distance2DSquared(const PointXY &p1, const PointXY &p2)
		{
			return (p1.x - p2.x)*(p1.x - p2.x) + (p1.y - p2.y)*(p1.y - p2.y);
		}
		/**\brief calculate distance between two 2D planes
		  */
		template<class PointXY>
		double distance2D(const PointXY &p1, const PointXY &p2)
		{
			return sqrt((p1.x - p2.x)*(p1.x - p2.x) + (p1.y - p2.y)*(p1.y - p2.y));
		}
		/**\brief insert sort
		  */
		template<class Iterator,class BinaryComp>
		void insert_sort(Iterator begin, Iterator end, BinaryComp pred) 
		{
			for (auto it = begin; it + 1 != end; ++it) {
				for (auto it2 = it + 1; it2 != begin; --it2) {
					if (pred(*it2, *(it2 - 1))) {
						auto temp = *(it2 - 1);
						*(it2 - 1) = *it2;
						*it2 = temp;
					}
					else {
						break;
					}
				}
			}
		}
		/**\brief merge sort, with nlogn temporal complexity
		*/
		template<class DataType,class BinaryComp>
		void merge_sort(std::vector<DataType> &list, BinaryComp pred)
		{
			if (list.size() <= 1)
				return;
			int half = list.size() / 2;
			std::vector<DataType> pre;
			std::vector<DataType> post;
			for (int i = 0; i < half; ++i) {
				pre.push_back(list[i]);
			}
			for (int i = half; i < list.size(); ++i) {
				post.push_back(list[i]);
			}
			merge_sort(pre, pred);
			merge_sort(post, pred);
			//�鲢
			auto it1 = pre.begin();
			auto it2 = post.begin();
			for (int i = 0; i < list.size(); ++i) {
				if (it1 == pre.end()) {
					list[i] = *it2++;
				}
				else if (it2 == post.end()) {
					list[i] = *it1++;
				}
				else {
					list[i] = pred(*it1, *it2) ? *it1++ : *it2++;
				}
			}
		}
	}
#endif
	/**\brief encapsulate numpy equivalent functionalities
	 * \note: this code has not been vectorized(it's implemented in loops) and thus not particularly fast
	  */
	namespace numpy
	{
		/**\brief limited version of broadcast like that in numpy-broadcast one number to a vector
		  *\param[in] val the number
		  *\param[in] num number of elements in the output vector
		  *\param[out] vec the resulting vector
		  */
		template<class NumType, class allocator = std::allocator<NumType> >
		void broadcast(const NumType &val, int num, std::vector<NumType, allocator> &vec)
		{
			vec.clear();
			vec.resize(num);
			for (auto &v : vec) { v = val; }
		}
		/**\brief limited version of broadcast like that in numpy-broadcast one vector to a 2d vector(matrix)
		  *\param[in] vec the input vector
		  *\param[in] repeats number of repeats in the output 2d vector
		  *\param[out] mat the resulting matrix
		  *\param[in] prepend true if broadcast along first dimension(namely repeat rows),false if broadcast along second dimesion(namely repeat columns)
		  */
		template<class NumType, class allocator = std::allocator<NumType> >
		void broadcast(const std::vector<NumType, allocator> &vec, int repeats, std::vector<std::vector<NumType, allocator> > &mat, bool prepend = true)
		{
			mat.clear();
			mat.resize(prepend ? repeats : vec.size());
			if (prepend) {
				for (int i = 0; i < repeats; ++i) {
					for (auto it = vec.begin(); it != vec.end(); ++it) {
						mat[i].push_back(*it);
					}
				}
			}
			else {
				for (int i = 0; i < mat.size(); ++i) {
					for (int j = 0; j < repeats; ++j) {
						mat[i].push_back(vec[i]);
					}
				}
			}
		}
		/**\brief limited version of broadcast like that in numpy-broadcast one number to a matrix
		  *\param[in] val the number
		  *\param[in] row,col the number of rows and cols in the resulting matrix
		  *\param[out] mat the resulting matrix
		  */
		template<class NumType, class allocator = std::allocator<NumType> >
		void broadcast(const NumType &val, int row, int col, std::vector<std::vector<NumType, allocator> > &mat)
		{
			mat.clear();
			mat.resize(row);
			for (int i = 0; i < row; ++i) {
				for (int j = 0; j < col; ++j) {
					mat[i].push_back(val);
				}
			}
		}
		/**\brief transpose the matrix
		  *\param[in] mat the input matrix
		  *\param[out] out_mat the transposed matrix
		  */
		template<class NumType, class allocator = std::allocator<NumType> >
		void transpose(const std::vector<std::vector<NumType, allocator> > &mat, std::vector<std::vector<NumType, allocator> > &out_mat)
		{
			out_mat.clear();
			if (mat.size() == 0) {
				return;
			}
			out_mat.resize(mat[0].size());
			for (int i = 0; i < mat[0].size(); ++i) {
				for (int j = 0; j < mat.size(); ++j) {
					out_mat[i].push_back(mat[j][i]);
				}
			}
		}
		template<class NumType,class allocator=std::allocator<NumType> >
		void range(std::vector<NumType, allocator> &result,NumType start, NumType end, NumType step = 1)
		{
			if (end == start || step == 0) {
				return;
			}
			if ((end - start)*step < 0) {//step is opposite to start-to-end direction
				throw std::invalid_argument("error from lyxutils::numpy::range()=>step=="+str_utils::type2str(step)+" won't take you from "+str_utils::type2str(start)+" to "+str_utils::type2str(end));
			}
			result.clear();
			for (NumType i = start; (i - end)*(i-end-step)>0; i += step) {
				result.push_back(i);
			}
		}
		template<class NumType, class allocator = std::allocator<NumType> >
		void range(std::vector<NumType, allocator> &result,NumType end)
		{
			range<NumType,allocator>(result, 0, end);
		}
		template<class NumType, class allocator = std::allocator<NumType> >
		void linspace(std::vector<NumType, allocator> &result,NumType start, NumType end, int num)
		{
			if (num < 1) {
				throw std::invalid_argument("illegal parameter num,num should be greater than 1(or equal to 1,only when start==end)");
			}
			if (num == 1) {
				if (start == end) {
					result.push_back(start);
					return;
				}
				else {
					throw std::invalid_argument("illegal parameter num,num should be greater than 1(or equal to 1,only when start==end)");
				}
			}
			NumType step = (end - start) / (num - 1);
			range(result,start, end, step);
			result.push_back(end);
		}
		/**\brief reshape the one dimensional vector into two dimensional
		  *\note: row times col should be equal to size of vec;if not the col will be automatically adjusted to satisfy the demand when row!=0 and size of vec is divisible by row
		  *\param[in] vec the 1d vector
		  *\param[in] row,col the output vector number of rows and cols
		  *\param[out] out_vec the resulting 2d vector
		  */
		template<class NumType, class allocator = std::allocator<NumType> >
		void reshape(const std::vector<NumType, allocator> &vec, int row, int col, std::vector<std::vector<NumType, allocator> > &out_vec)
		{
			if (row*col != vec.size()) {
				if (row != 0 && vec.size() % row == 0 && col==-1) {
					col = vec.size() / row;
				}
				else if(col !=0 && vec.size() % col ==0 && row==-1){
				    row = vec.size() / col;
				}
				else {
					throw std::invalid_argument("row times col is not equal to size of vec!");
				}
			}
			out_vec.clear();
			out_vec.resize(row);
			for (int i = 0; i < row; ++i) {
				for (int j = 0; j < col; ++j) {
					out_vec[i].push_back(vec[i*col + j]);
				}
			}
		}
		/**\brief expand the two dimensional vector into 1d one
		  *\param[in] vec the input two d vector
		  *\param[out] out_vec the expanded vector
		  */
		template<class NumType, class allocator = std::allocator<NumType> >
		void flat(const std::vector<std::vector<NumType, allocator> > &vec, std::vector<NumType, allocator> &out_vec)
		{
			if (vec.size() == 0 || vec[0].size() == 0) {
				return;
			}
			out_vec.clear();
			for (int i = 0; i < vec.size(); ++i) {
				for (int j = 0; j < vec[i].size(); ++j) {
					out_vec.push_back(vec[i][j]);
				}
			}
		}
		template<class NumType, class allocator = std::allocator<NumType> >
		void print(const std::vector<NumType, allocator> &vec, bool column = false, bool keepSameLine = false, bool verbose = true)
		{
			std::cout << "[";
			if (vec.size() > 0)
			{
				std::cout << vec[0];
				if (!column) {
					for (auto it = vec.begin() + 1; it != vec.end(); ++it) {
						std::cout << ",\t" << *it;
					}
				}
				else {
					for (auto it = vec.begin() + 1; it != vec.end(); ++it) {
						std::cout << "," << std::endl << " " << *it;
					}
				}
			}
			std::cout << "]";
			if (verbose) {
				std::cout << ",number of elements:" << vec.size() << ",type of element:" << std::string(typeid(NumType).name());
			}
			if (!keepSameLine) {
				std::cout << std::endl;
			}
		}
		template<class NumType, class allocator = std::allocator<NumType> >
		void print(const std::vector<std::vector<NumType, allocator> > &mat, bool keepSameLine = false, bool verbose = true)
		{
			std::cout << "[";
			for (auto it = mat.begin(); it != mat.end(); ++it) {
				if (it != mat.begin()) {
					std::cout << " ";
				}
				print(*it, false, it == (mat.end() - 1),false);
			}
			std::cout << "]";
			if (verbose) {
				std::cout << ",shape:" << mat.size() << "x" << mat[0].size() << ",type of element:" << std::string(typeid(NumType).name());
			}
			if (!keepSameLine) {
				std::cout << std::endl;
			}
		}
		/**\brief get the summation
		  */
		template<class NumType,class allocator=std::allocator<NumType> >
		NumType sum(const std::vector<NumType,allocator> &vec)
		{
			NumType result = {0};
			for (auto it = vec.begin(); it != vec.end(); ++it)result += (*it);
			return result;
		}
		/**\brief get the summation along one axis
		  *\param[in] vec input 2D qrray(first dimension is row)
		  *\param[out] out_vec output 1D array-the summation
		  *\param[in] axis the direction of summation(default -1,namely add all elements up)
		  *\note:axis=0->sum along row direction; axis=1->sum along column direction; axis=-1->sum all the elements
		  *\example:
		  *\std::vector<float> sum;
		  *\lyxutils::stats::sum(input_vec,sum,axis=1);//sum all the rows
		*/
		template<class NumType, class allocator = std::allocator<NumType> >
		void sum(const std::vector<std::vector<NumType, allocator> > &mat, std::vector<NumType, allocator> &out_vec, int axis = -1)
		{
			if (mat.size() == 0 || mat[0].size() == 0) {
				out_vec.clear();
				out_vec.push_back({0});
				return;
			}
			if (axis == 0 || axis == -1)
			{
				out_vec.resize(mat[0].size());
				for (NumType &s : out_vec)s = {0};
				for (int i = 0; i < out_vec.size(); ++i)
				{
					for (auto &row : mat) {
						out_vec[i] += row[i];
					}
				}
				if (axis == -1) {
					NumType result = sum(out_vec);
					out_vec.clear();
					out_vec.push_back(result);
				}
			}
			else if (axis == 1)
			{
				out_vec.resize(mat.size());
				for (NumType &s : out_vec)s = {0};
				for (int i = 0; i < out_vec.size(); ++i)
				{
					for (auto &s : mat[i]) {
						out_vec[i] += s;
					}
				}
			}
			else {
				throw std::invalid_argument("invalid axis parameter!");
			}
		}
		/**\brief get average of 1D vector
		  */
		template<class NumType, class allocator = std::allocator<NumType> >
		NumType mean(const std::vector<NumType, allocator> &vec)
		{
			return sum(vec) / vec.size();
		}
		/**\brief get averages of 2D array, refer to lyxutils::stats::sum
		  */
		template<class NumType, class allocator = std::allocator<NumType> >
		void mean(const std::vector<std::vector<NumType, allocator> > &mat, std::vector<NumType, allocator> &out_vec, int axis = -1)
		{
			sum(mat, out_vec, axis);
			if (axis == -1) {
				out_vec[0] /= (mat.size()*mat[0].size());
			}
			else if (axis == 0) {
				for (int i = 0; i < out_vec.size(); ++i) {
					out_vec[i] /= mat.size();
				}
			}
			else if (axis == 1) {
				for (int i = 0; i < out_vec.size(); ++i) {
					out_vec[i] /= mat[i].size();
				}
			}
			else {
				throw std::invalid_argument("invalid axis parameter!");
			}
		}
		/**\brief get the standard deviation of 1D vector
		  *\param[in] vec input vector
		  *\param[in] ddof degree of freedom
		  *\std=sqrt(sum((xi-xbar)^2)/(n-ddof)),where xi represents each ith element, xbar is the average x value
		  *\n is number of elements and ddof is degree of freedom
		  */
		template<class NumType, class allocator = std::allocator<NumType> >
		NumType std(const std::vector<NumType, allocator> &vec, int ddof = 1)
		{
			NumType xbar = mean(vec);
			int n = vec.size();
			std::vector<NumType, allocator> xi2;
			apply(vec, xi2, [ddof, xbar, n](NumType x)->float {return (x - xbar)*(x - xbar) / (n - ddof); });
			return sqrt(sum(xi2));
		}
		/**\brief get the standard deviation of 2D matrix
		  *\param[in] mat input matrix
		  *\param[out] out_vec output
		  *\param[in] ddof degree of freedom
		  *\std=sqrt(sum((xi-xbar)^2)/(n-ddof)),where xi represents each ith element, xbar is the average x value
		  *\n is number of elements and ddof is degree of freedom
		  */
		template<class NumType, class allocator = std::allocator<NumType> >
		void std(const std::vector<std::vector<NumType, allocator> > &mat, std::vector<NumType, allocator> &out_vec, int ddof = 1, int axis = -1)
		{
			if (mat.size() == 0 || mat[0].size() == 0) {
				out_vec.clear();
				out_vec.push_back(0);
				return;
			}
			if (axis == -1) {
				std::vector<NumType, allocator> xbar;
				mean(mat, xbar);
				NumType m = xbar[0];
				int n = mat.size()*mat[0].size();
				NumType sqrSum = 0;
				for (int i = 0; i < mat.size(); ++i) {
					for (int j = 0; j < mat[i].size(); ++j) {
						sqrSum += (mat[i][j] - m)*(mat[i][j] - m) / (n - ddof);
					}
				}
				out_vec.clear();
				out_vec.push_back(sqrt(sqrSum));
			}
			else if (axis == 0) {
				out_vec.resize(mat[0].size());
				std::vector<NumType, allocator> rowMean;
				mean(mat, rowMean, 0);
				int n = mat.size();
				for (auto &s : out_vec) { s = 0; }
				for (int i = 0; i < out_vec.size(); ++i) {
					for (int j = 0; j < mat.size();++j) {
						out_vec[i] += (mat[j][i] - rowMean[i])*(mat[j][i] - rowMean[i]) / (n - ddof);
					}
				}
				for (auto &s : out_vec) { s = sqrt(s); }
			}
			else if (axis == 1) {
				out_vec.resize(mat.size());
				std::vector<NumType, allocator> colMean;
				mean(mat, colMean, 1);
				int n = mat[0].size();
				for (auto &s : out_vec) { s = 0; }
				for (int i = 0; i < out_vec.size(); ++i) {
					for (int j = 0; j < mat[0].size(); ++j) {
						out_vec[i] += (mat[i][j] - colMean[i])*(mat[i][j] - colMean[i]) / (n - ddof);
					}
				}
				for (auto &s : out_vec) { s = sqrt(s); }
			}
			else {
				throw std::invalid_argument("invalid axis parameter!");
			}
		}
		/**\brief add two vectors
		  *\param[in] vec1,vec2 the input vector of the same size
		  *\param[out] out_vec sum of input vectors
		  */
		template<class NumType, class allocator = std::allocator<NumType> >
		void add(const std::vector<NumType, allocator> &vec1, const std::vector<NumType, allocator> &vec2, std::vector<NumType, allocator> &out_vec)
		{
			if (vec1.size() == 0 && vec2.size() == 0) {
                out_vec.clear();
			    //do nothing if sizes of both vectors are zero
			}
			else if (vec1.size() == vec2.size()) {
                out_vec.clear();
				for (auto it1 = vec1.begin(), it2 = vec2.begin(); it1 != vec1.end(); ++it1, ++it2) {
					out_vec.push_back((*it1)+(*it2));
				}
			}
			else {
				throw std::runtime_error("size of vec1 is not compatible with that of vec2!");
			}
		}
		/**\brief add two matrices
		  *\param[in] matrix1,matrix2 the input matrices of the same size
		  *\param[out] out_mat sum of input matrices
		  */
		template<class NumType, class allocator = std::allocator<NumType> >
		void add(const std::vector<std::vector<NumType, allocator> > &matrix1, const std::vector<std::vector<NumType, allocator> > &matrix2, std::vector<std::vector<NumType, allocator> > &out_mat)
		{
			if (matrix1.size() == 0 && matrix2.size() == 0) {
                out_mat.clear();
			}
			else if (matrix1.size() == matrix2.size()) {
                out_mat.clear();
				for (auto it1 = matrix1.begin(), it2 = matrix2.begin(); it1 != matrix1.end(); ++it1, ++it2) {
					out_mat.push_back(std::vector<NumType, allocator>());
					add(*it1, *it2, *(out_mat.end() - 1));
				}
			}
			else {
				throw std::runtime_error("size of matrix1 is not compatible with that of matrix2!");
			}
		}
		/**\brief minus two vectors
		  *\param[in] vec1,vec2 the input vector of the same size
		  *\param[out] out_vec result of input vectors
		  */
		template<class NumType, class allocator = std::allocator<NumType> >
		void minus(const std::vector<NumType, allocator> &vec1, const std::vector<NumType, allocator> &vec2, std::vector<NumType, allocator> &out_vec)
		{
			if (vec1.size() == 0 && vec2.size() == 0) {
				//do nothing if sizes of both vectors are zero
                out_vec.clear();
			}
			else if (vec1.size() == vec2.size()) {
                out_vec.clear();
				for (auto it1 = vec1.begin(), it2 = vec2.begin(); it1 != vec1.end(); ++it1, ++it2) {
					out_vec.push_back((*it1) - (*it2));
				}
			}
			else {
				throw std::runtime_error("size of vec1 is not compatible with that of vec2!");
			}
		}
		/**\brief minus two matrices
		  *\param[in] matrix1,matrix2 the input matrices of the same size
		  *\param[out] out_mat result of input matrices
		  */
		template<class NumType, class allocator = std::allocator<NumType> >
		void minus(const std::vector<std::vector<NumType, allocator> > &matrix1, const std::vector<std::vector<NumType, allocator> > &matrix2, std::vector<std::vector<NumType, allocator> > &out_mat)
		{
			if (matrix1.size() == 0 && matrix2.size() == 0) {
                out_mat.clear();
			}
			else if (matrix1.size() == matrix2.size()) {
                out_mat.clear();
				for (auto it1 = matrix1.begin(), it2 = matrix2.begin(); it1 != matrix1.end(); ++it1, ++it2) {
					out_mat.push_back(std::vector<NumType, allocator>());
					minus(*it1, *it2, *(out_mat.end() - 1));
				}
			}
			else {
				throw std::runtime_error("size of matrix1 is not compatible with that of matrix2!");
			}
		}
		/**\brief compute the dot product of two vectors
		  *\param[in] vec1,vec2 the input vector of the same size
		  *\param[out] out_vec the output vector, which is the dot product of input vectors
		  */
		template<class NumType,class allocator=std::allocator<NumType> >
		void dot(const std::vector<NumType, allocator> &vec1, const std::vector<NumType, allocator> &vec2, std::vector<NumType, allocator> &out_vec)
		{
			if (vec1.size() == 0 && vec2.size() == 0) {
				//do nothing if sizes of both vectors are zero
                out_vec.clear();
			}
			else if (vec1.size() == vec2.size()) {
                out_vec.clear();
				for (auto it1 = vec1.begin(), it2 = vec2.begin(); it1 != vec1.end(); ++it1, ++it2) {
					out_vec.push_back((*it1)*(*it2));
				}
			}
			else {
				throw std::runtime_error("size of vec1 is not compatible with that of vec2!");
			}
		}
		/**\brief compute the dot product of two matrices
		*\param[in] matrix1,matrix2 the input matrices of the same size
		*\param[out] out_mat the output matrix, which is the dot product of input matrices
		*/
		template<class NumType, class allocator = std::allocator<NumType> >
		void dot(const std::vector<std::vector<NumType, allocator> > &matrix1, const std::vector<std::vector<NumType, allocator> > &matrix2, std::vector<std::vector<NumType, allocator> > &out_mat)
		{
			if (matrix1.size() == 0 && matrix2.size() == 0) {
                out_mat.clear();
			}
			else if (matrix1.size() == matrix2.size()) {
                out_mat.clear();
				for (auto it1 = matrix1.begin(), it2 = matrix2.begin(); it1 != matrix1.end(); ++it1, ++it2) {
					out_mat.push_back(std::vector<NumType, allocator>());
					dot(*it1, *it2, *(out_mat.end() - 1));
				}
			}
			else {
				throw std::runtime_error("size of matrix1 is not compatible with that of matrix2!");
			}
		}
		/**\brief compute the inner product of two vectors
		  *\param[in] vec1,vec2 the input vector of the same size
		  */
		template<class NumType, class allocator = std::allocator<NumType> >
		NumType multiply(const std::vector<NumType, allocator> &vec1, const std::vector<NumType, allocator> &vec2)
		{
			NumType sum = {0};
			if (vec1.size() == 0 || vec2.size() == 0) {
				//do nothing if sizes of both vectors are zero
			}
			else if (vec1.size() == vec2.size()) {
				for (auto it1 = vec1.begin(), it2 = vec2.begin(); it1 != vec1.end(); ++it1, ++it2) {
					sum += (*it1)*(*it2);
				}
			}
			else {
				throw std::runtime_error("size of vec1 is not compatible with that of vec2!");
			}
			return sum;
		}
		/**\brief compute the matrix product of two vectors, consider first vector as column vector matrix and second as row vector matrix
		  *\param[in] vec1,vec2 the input vector(size of vec1=m,size of vec2=n)
		  *\param[out] result the resulting matrix of the matrix product(size:m*n)
		  */
		template<class NumType, class allocator = std::allocator<NumType> >
		void multiply(const std::vector<NumType, allocator> &vec1, const std::vector<NumType, allocator> &vec2, std::vector<std::vector<NumType, allocator> > &result)
		{
			result.clear();
			result.resize(vec1.size());
			for (int i = 0; i < result.size(); ++i) {
				result[i].resize(vec2.size());
				for (int j = 0; j < vec2.size(); ++j) {
					result[i][j] = vec1[i] * vec2[j];
				}
			}
		}
		/**\brief compute the matrix product of two matrices(underlying assumption:size of all second dimensions are the same)
		  *\param[in] matrix1,matrix2 the input matrices(cols of matrix1 should equal rows of matrix2)
		  *param[out] out_mat the resultant matrix product of the two input matrices
		  */
		template<class NumType, class allocator = std::allocator<NumType> >
		void multiply(const std::vector<std::vector<NumType, allocator> > &matrix1, const std::vector<std::vector<NumType, allocator> > &matrix2, std::vector<std::vector<NumType, allocator> > &out_mat)
		{
		    if(matrix1.size()==0||matrix1[0].size()==0||
		       matrix2.size()==0||matrix2[0].size()==0){
		        throw std::runtime_error("size of matrix1 and matrix2 cannot be zero!");
		    }
		    if(matrix1[0].size()!=matrix2.size()){
		        throw std::logic_error("columns in matrix1 should equal rows in matrix2 when matrix1 x matrix2!");
		    }
			out_mat.clear();
			for (int row = 0; row < matrix1.size(); ++row) {
				out_mat.push_back(std::vector<NumType,allocator>());
				for (int col = 0; col < matrix2[0].size(); ++col) {
					out_mat[row].push_back(0);
					NumType &sum = out_mat[row][col];
					auto it1 = matrix1[row].begin();
					auto it2 = matrix2.begin();
					for (; it1 != matrix1[row].end(); ++it1, ++it2) {
						sum += *it1*(*it2)[col];
					}
				}
			}
		}
		/**\brief apply a function to each element of a vector,eg.dst=f(src)
		  *\param[in] src the input vector
		  *\param[out] dst the output vector
		  *\param,[in] f the applied function,f should take in a NumType1 argument and return a NumType2 value
		  */
		template<class NumType1, class NumType2, class allocator1 = std::allocator<NumType1>,class allocator2=std::allocator<NumType2>, class Func>
		void apply(const std::vector<NumType1, allocator1> &src, std::vector<NumType2, allocator2> &dst,Func f)
		{
			dst.clear();
			dst.resize(src.size());
			auto it1 = src.begin();
			auto it2 = dst.begin();
			for (; it1 != src.end(); ++it1, ++it2) {
				*it2 = f(*it1);
			}
		}
		/**\brief apply a function to each element of a vector,eg.f(src)
		  *\param[in] src the input vector
		  *\param,[in] f the applied function,f should take in a NumType argument
		  */
		template<class NumType, class allocator = std::allocator<NumType>, class Func>
		void apply(std::vector<NumType, allocator> &src, Func f)
		{
			auto it = src.begin();
			for (; it != src.end(); ++it) {
				f(*it);
			}
		}
		/**\brief apply a function to each element of a matrix,eg.dst=f(src)
		  *\param[in] src the input matrix
		  *\param[out] dst the output matrix
		  *\param,[in] f the applied function,f should take in a NumType1 argument and return a NumType2 value
		  */
		template<class NumType1, class NumType2, class allocator1 = std::allocator<NumType1>, class allocator2 = std::allocator<NumType2>, class Func>
		void apply(const std::vector<std::vector<NumType1, allocator1> > &src, std::vector<std::vector<NumType2, allocator2> > &dst, Func f)
		{
		    if(src.size()==0||src[0].size()==0){
		        return;
		    }
			dst.clear();
			dst.resize(src.size());
			auto it1 = src.begin();
			auto it2 = dst.begin();
			for (; it1 != src.end(); ++it1, ++it2) {
				apply(*it1, *it2, f);
			}
		}
		/**\brief apply a function to each element of a matrix,eg.f(src)
		  *\param[in] src the input matrix
		  *\param,[in] f the applied function,f should take in a NumType argument
		  */
		template<class NumType, class allocator = std::allocator<NumType>, class Func>
		void apply(std::vector<std::vector<NumType, allocator> > &src, Func f)
		{
			auto it = src.begin();
			for (; it != src.end(); ++it) {
				apply(*it, f);
			}
		}
		/**\brief find the elements in a vector that meet a certain requirement
		  *\param[in] src the given vector
		  *\param[out] indices the indices of elements that satisfy the requirement(s)
		  *\param[in] f the condition function, this function takes in a NumType argument and returns a boolean value that decides whether the NumType element satisfy the requirement(s)
		  */
		template<class NumType, class allocator = std::allocator<NumType>, class Func>
		void where(const std::vector<NumType, allocator> &src, std::vector<int> &indices, Func f)
		{
			indices.clear();
			for (int i = 0; i < src.size(); ++i) {
				if (f(src[i])) {
					indices.push_back(i);
				}
			}
		}
		/**\brief find the elements in a matrix that meet a certain requirement
		  *\param[in] src the given matrix
		  *\param[out] indix_x,index_y the indices of elements that satisfy the requirement(s)
		  *\param[in] f the condition function, this function takes in a NumType argument and returns a boolean value that decides whether the NumType element satisfy the requirement(s)
		  */
		template<class NumType, class allocator = std::allocator<NumType>, class Func>
		void where(const std::vector<std::vector<NumType, allocator> > &src, std::vector<int> &index_x, std::vector<int> &index_y, Func f)
		{
			index_x.clear();
			index_y.clear();
			for (int i = 0; i < src.size(); ++i) {
				for (int j = 0; j < src[i].size(); ++j) {
					if (f(src[i][j])) {
						index_x.push_back(i);
						index_y.push_back(j);
					}
				}
			}
		}
		/**\brief find the elements in a vector that meet a certain requirement and apply function to it
		  *\param[in] src the given vector
		  *\param[in] f the condition function, this function takes in a NumType argument and returns a boolean value that decides whether the NumType element satisfy the requirement(s)
		  *\param[in] pro the applied function
		  */
		template<class NumType, class allocator = std::allocator<NumType>, class Pred,class Process>
		void apply_where(std::vector<NumType, allocator> &src, Pred f, Process f2)
		{
			for (int i = 0; i < src.size(); ++i) {
				if (f(src[i])) {
					f2(src[i]);
				}
			}
		}
		/**\brief find the elements in a matrix that meet a certain requirement and apply function to it
		  *\param[in] src the given matrix
		  *\param[in] f the condition function, this function takes in a NumType argument and returns a boolean value that decides whether the NumType element satisfy the requirement(s)
		  *\param[in] pro the applied function
		  */
		template<class NumType, class allocator = std::allocator<NumType>, class Pred, class Process>
		void apply_where(std::vector<std::vector<NumType, allocator> > &src, Pred f, Process f2)
		{
			for (int i = 0; i < src.size(); ++i) {
				for (int j = 0; j < src[i].size(); ++j) {
					if (f(src[i][j])) {
						f2(src[i][j]);
					}
				}
			}
		}
        /**
         * get the covariance matrix of a matrix
         * @tparam NumType type of element
         * @tparam allocator allocator of vector
         * @param mat input matrix m*n=>represents m features each of n dimensions
         * @param covariance output covariance matrix n*n
         */
        template<class NumType, class allocator = std::allocator<NumType> >
        void cov(const std::vector<std::vector<NumType, allocator> > &mat, std::vector<std::vector<NumType, allocator> > &covariance){
            std::vector<NumType, allocator> m;
            mean(mat,m,0);
            std::vector<std::vector<NumType,allocator> > m_ext;
            broadcast(m,mat.size(),m_ext);
            for(int i=0;i<mat.size();++i){
                for(int j=0;j<mat[i].size();++j){
                    m_ext[i][j]=mat[i][j]-m_ext[i][j];
                }
            }
            std::vector<std::vector<NumType,allocator> > m_trans;
            transpose(m_ext,m_trans);
            multiply(m_trans,m_ext,covariance);
        }
	}
	/**\brief encapsulate some os related functionalities
	 */
	namespace os{
	    //the maximum path length returned by pwd()
	    const int __MAX_PATH_SIZE=256;
	    /**\brief check if a path is an empty string and throw exception for this reason
	     */
	    inline void __check_path_validity(std::string path){
            if(path=="")throw std::logic_error("logic error from os::__check_path_validity:path should not be an empty string");
	    }
	    /**\check if file name use relative path or absolute path
	     */
	    bool __isrelative(std::string file);
	    /**\brief check if fileName exists
	     */
	    bool exist(std::string fileName);
	    /**\brief check if fileName is a regular file
	     */
	    bool file(std::string fileName);
	    /**\brief check if path exists and is a directory
	     */
	    bool directory(std::string path);
	    /**\brief present working directory
	     */
	    std::string pwd();
	    /**\brief get parent path of a given path/file
	     */
        std::string parent(std::string path);
        /**\brief list all children of directory path. if path does not exist or is not a directory return empty list
         * this function throws std::runtime_error when open folder failed
         * param[in] path the parent folder to be listed;
         * param[in] flag 0:list all files,1:only regular files,2:only directories
         * return children's absolute name(with path and suffix)
         */
        std::vector<std::string> ls(std::string path=pwd(), int flag=0);
        /**ls all directories contained in path.if path does not exist or is not a directory return empty list
         * return children directories' absolute name(with path and suffix)
         */
        std::vector<std::string> lsdir(std::string path=pwd());
        /**\brief list all regular files contained in path.if path does not exist or is not a directory return empty list
         * return children files' absolute name(with path and suffix)
         */
        std::vector<std::string> lsfile(std::string path=pwd());
        /**\brief find the reletive path of target path(directory or file) from base path(directory or file)
         */
        std::string relpath(std::string target,std::string base=pwd());
        /**\brief get the absolute path name of fileName
         */
        std::string abspath(std::string fileName);
        /**\brief get the file name without path(remove everything before '/' or '\\')
         * this function does not check the existence of given file or whether it is a regular file
         */
        std::string fullname(std::string fileName);
        /**\brief get the base name of an existing file without path and suffix
         * consider the string after last dot(.) as suffix
         * this function does not check the existence of given file or whether it is a regular file
         */
        std::string basename(std::string fileName);
        /**\brief get the suffix of a file.
         * this function does not check the existence of given file or whether it is a regular file
         */
        std::string suffix(std::string fileName);
        /**\brief change current directory.support relative path
         */
        bool cd(std::string target);
        /**\brief make a directory.
         * this function doesn't set permissions under windows and set permissions=775 in linux
         * if this is not what you want,use function ::mkdir from stat.h or use system command
         * function touch,rmdir and rm are deleted because create a file is easy using ofstream;
         * remove files and directories are dangerous operations that ain't as frequently used in my programming
         * use system("touch file"),system("rmdir dir") and system("rm file") instead
         */
        bool mkdir(std::string path, bool recursive=false);
	}
	/**\brief encapsulate some useful geometry-related data structure and functions
	 */
	namespace geometry{
	    /**\brief get the return value type of the result of a binary operator
	     * @tparam A type of the first operand of the binary operator
	     * @tparam B type of the second operand of the binary operator
	     */
	    template<typename A,typename B>
	    struct __BiOperatorType{
	    private:
	        static A __A();
	        static B __B();
	    public:
	        typedef decltype(__A()+__B()) AddType;
            typedef decltype(__A()-__B()) MinusType;
            typedef decltype(__A()*__B()) MultiType;
            typedef decltype(__A()/__B()) DivideType;
	    };
        //forward declaration
        template<typename T1,typename T2>
        struct __RTDAdd;
        template<typename T1,typename T2>
        struct __RTDMinus;
        template<typename T1,typename T2>
        struct __RTDMultiply;
        template<typename T1,typename T2>
        struct __RTDDivide;
        /**
         * convenient 3d vector for any value type(int,float,double...)
         */
        template<typename T>
        class Vector3{
        public:
            T x,y,z;
            typedef T value_type;
            Vector3():Vector3({}){}//initialized by zeros
            template<typename T1>
            Vector3(const Vector3<T1> &other):x(other.x),y(other.y),z(other.z){}
            Vector3(const T &_x, const T &_y, const T &_z):x(_x),y(_y),z(_z){}
            Vector3(std::initializer_list<T> init_list){
                auto it=init_list.begin();
                for(int i=0;i<3;++i){
                    if(it!=init_list.end()){
                        (*this)[i]=*it;
                        ++it;
                    }
                    else{
                        (*this)[i]=0;
                    }
                }
            }
            Vector3(const std::vector<T> &v, int start=0){
                auto it=v.begin()+start;
                for(int i=0;i<3;++i){
                    if(start+i<v.size()){
                        (*this)[i]=*it;
                        ++it;
                    }
                    else{
                        (*this)[i]=0;
                    }
                }
            }
            template<typename T1>
            Vector3<T> operator=(const Vector3<T1> &v){
                this->x=v.x;
                this->y=v.y;
                this->z=v.z;
                return *this;
            }
            T& operator[](int index){
                if(index==0)return x;
                if(index==1)return y;
                if(index==2)return z;
                throw std::out_of_range("Vector3 index out of range!");
            }
            //as the standard practice in stl did, add a const_reference access operator[]
            const T& operator[](int index)const{
                if(index==0)return x;
                if(index==1)return y;
                if(index==2)return z;
                throw std::out_of_range("Vector3 index out of range!");
            }
            Vector3<T> operator-(){
                return {-x,-y,-z};
            }
            template<typename T1>
            auto operator+(const T1 &t)-> typename __RTDAdd<Vector3<T>,T1>::ret_type{
                return __RTDAdd<Vector3<T>,T1>().add(*this,t);
            }
            template<typename T1>
            auto operator-(const T1 &t)-> typename __RTDMinus<Vector3<T>,T1>::ret_type{
                return __RTDMinus<Vector3<T>,T1>().minus(*this,t);
            }
            template<typename T1>
            auto operator*(const T1 &t)-> typename __RTDMultiply<Vector3<T>,T1>::ret_type{
                return __RTDMultiply<Vector3<T>,T1>().multiply(*this,t);
            }
            template<typename T1>
            auto operator/(const T1 &t)-> typename __RTDDivide<Vector3<T>,T1>::ret_type{
                return __RTDDivide<Vector3<T>,T1>().divide(*this,t);
            }
            //friend operator +,-,*,/ take parameters T1 and Vector3<T2> not T1 and T2
            // in case capture all +,-,*,/ operations
            template<typename T1,typename T2>
            friend auto operator+(const T1 &t1,const Vector3<T2> &t2)->typename __RTDAdd<T1,Vector3<T2> >::ret_type;
            template<typename T1,typename T2>
            friend auto operator-(const T1 &t1,const Vector3<T2> &t2)->typename __RTDMinus<T1,Vector3<T2> >::ret_type;
            template<typename T1,typename T2>
            friend auto operator*(const T1 &t1,const Vector3<T2> &t2)->typename __RTDMultiply<T1,Vector3<T2> >::ret_type;
            template<typename T1,typename T2>
            friend auto operator/(const T1 &t1,const Vector3<T2> &t2)->typename __RTDDivide<T1,Vector3<T2> >::ret_type;
            template<typename T1>
            void operator +=(const T1 right){
                auto res=__RTDAdd<Vector3<T>,T1>().add(*this,right);
                x=res.x;
                y=res.y;
                z=res.z;
            }
            template<typename T1>
            void operator -=(const T1 right){
                auto res=__RTDMinus<Vector3<T>,T1>().minus(*this,right);
                x=res.x;
                y=res.y;
                z=res.z;
            }
            template<typename T1>
            void operator *=(const T1 right){
                auto res=__RTDMultiply<Vector3<T>,T1>().multiply(*this,right);
                x=res.x;
                y=res.y;
                z=res.z;
            }
            template<typename T1>
            void operator /=(const T1 right){
                auto res=__RTDDivide<Vector3<T>,T1>().divide(*this,right);
                x=res.x;
                y=res.y;
                z=res.z;
            }
            //dot product
            template<typename T1>
            auto dot(const Vector3<T1> &v)const -> decltype(x*v.x){
                return x*v.x+y*v.y+z*v.z;
            }
            //cross product
            template<typename T1>
            auto cross(const Vector3<T1> &v)const-> Vector3<decltype(x*v.y)>{
                typedef decltype(x*v.x) element_type;
                return Vector3<element_type>(y*v.z-z*v.y,z*v.x-x*v.z,x*v.y-y*v.x);
            }
            //mold,rank=0 means L0 norm=>L0(x)=||x||_0=#(i|x_i!=0)=>represents number of non-zeros in vector
            T mold(double rank=2)const{
                if(rank==0)return T((x!=0)+(y!=0)+(z!=0));
                return std::pow(std::pow(x,rank)+std::pow(y,rank)+std::pow(z,rank),1/rank);
            }
            template<typename T1>
            friend std::ostream& operator<<(std::ostream &out,const Vector3<T1> &v);
            /**\brief compare the equaty of two Vector3<T>
             * be careful about the tolerance caused by floating point number rounding error
             */
            template<typename T1>
            bool equal(const Vector3<T1> &t, typename __RTDMinus<T,T1>::ret_type tolerance={0}){
                //using abs(like following) will not produce correct result.don't know why
//                return abs(x-t.x)<=tolerance && abs(y-t.y)<=tolerance && abs(z-t.z)<=tolerance;
                Vector3<decltype(tolerance)> diff=*this-t;
                return diff.x<=tolerance && -diff.x<=tolerance &&
                       diff.y<=tolerance && -diff.y<=tolerance &&
                       diff.z<=tolerance && -diff.z<=tolerance;
            }
            template<typename T1>
            Vector3<T1> asType(){
                return Vector3<T1>(*this);
            }
            void fromStringVector(const std::vector<std::string> &v, int start= 0){
                auto it=v.begin()+start;
                for(int i=0;i<3;++i){
                    if(start+i<v.size()){
                        std::stringstream ss(*it);
                        ss>>(*this)[i];
                        ++it;
                    }
                    else{
                        (*this)[i]=0;
                    }
                }
            }
            std::vector<T> toStdVector(){
                return std::vector<T>({x,y,z});
            }
        };
        //helper structs of Vector3,only the partial specialization forms are used
        /**
         * __RTDAdd-return type deduced add operation.
         * specialize this template to define add operation for two customized type
         * @tparam T1 first operand type of add operation
         * @tparam T2 second operand type of add operation
         */
        template<typename T1,typename T2>
        struct __RTDAdd{
            typedef typename __BiOperatorType<T1,T2>::AddType ret_type;
            auto add(const T1 &t1,const T2 &t2)->ret_type{
                return t1+t2;
            }
        };
        template<typename T1,typename T2>
        struct __RTDAdd<Vector3<T1>,T2>{
            typedef Vector3<typename __BiOperatorType<T1,T2>::AddType> ret_type;
            auto add(const Vector3<T1> &t1,const T2 &t2)->ret_type{
                return ret_type(t1.x+t2,t1.y+t2,t1.z+t2);
            }
        };
        template<typename T1,typename T2>
        struct __RTDAdd<T1,Vector3<T2> >{
            typedef Vector3<typename __BiOperatorType<T1,T2>::AddType> ret_type;
            auto add(const T1 &t1,const Vector3<T2> &t2)->ret_type{
                return ret_type(t1+t2.x,t1+t2.y,t1+t2.z);
            }
        };
        template<typename T1,typename T2>
        struct __RTDAdd<Vector3<T1>,Vector3<T2> >{
            typedef Vector3<typename __BiOperatorType<T1,T2>::AddType> ret_type;
            auto add(const Vector3<T1> &t1,const Vector3<T2> &t2)->ret_type{
                return ret_type(t1.x+t2.x,t1.y+t2.y,t1.z+t2.z);
            }
        };
        /**
         * __RTDMinus-return type deduced minus operation.
         * specialize this template to define minus operation for two customized type
         * @tparam T1 first operand type of minus operation
         * @tparam T2 second operand type of minus operation
         */
        template<typename T1,typename T2>
        struct __RTDMinus{
            typedef typename __BiOperatorType<T1,T2>::MinusType ret_type;
            auto minus(const T1 &t1,const T2 &t2)->ret_type{
                return t1-t2;
            }
        };
        template<typename T1,typename T2>
        struct __RTDMinus<Vector3<T1>,T2>{
            typedef Vector3<typename __BiOperatorType<T1,T2>::MinusType> ret_type;
            auto minus(const Vector3<T1> &t1,const T2 &t2)->ret_type{
                return ret_type(t1.x-t2,t1.y-t2,t1.z-t2);
            }
        };
        template<typename T1,typename T2>
        struct __RTDMinus<T1,Vector3<T2> >{
            typedef Vector3<typename __BiOperatorType<T1,T2>::MinusType> ret_type;
            auto minus(const T1 &t1,const Vector3<T2> &t2)->ret_type{
                return ret_type(t1-t2.x,t1-t2.y,t1-t2.z);
            }
        };
        template<typename T1,typename T2>
        struct __RTDMinus<Vector3<T1>,Vector3<T2> >{
            typedef Vector3<typename __BiOperatorType<T1,T2>::MinusType> ret_type;
            auto minus(const Vector3<T1> &t1,const Vector3<T2> &t2)->ret_type{
                return ret_type(t1.x-t2.x,t1.y-t2.y,t1.z-t2.z);
            }
        };
        /**
         * __RTDMultiply-return type deduced multiply operation.
         * specialize this template to define multiply operation for two customized type
         * @tparam T1 first operand type of multiply operation
         * @tparam T2 second operand type of multiply operation
         */
        template<typename T1,typename T2>
        struct __RTDMultiply{
            typedef typename __BiOperatorType<T1,T2>::MultiType ret_type;
            auto multiply(const T1 &t1,const T2 &t2)->ret_type{
                return t1*t2;
            }
        };
        template<typename T1,typename T2>
        struct __RTDMultiply<Vector3<T1>,T2>{
            typedef Vector3<typename __BiOperatorType<T1,T2>::MultiType> ret_type;
            auto multiply(const Vector3<T1> &t1,const T2 &t2)->ret_type{
                return ret_type(t1.x*t2,t1.y*t2,t1.z*t2);
            }
        };
        template<typename T1,typename T2>
        struct __RTDMultiply<T1,Vector3<T2> >{
            typedef Vector3<typename __BiOperatorType<T1,T2>::MultiType> ret_type;
            auto multiply(const T1 &t1,const Vector3<T2> &t2)->ret_type{
                return ret_type(t1*t2.x,t1*t2.y,t1*t2.z);
            }
        };
        template<typename T1,typename T2>
        struct __RTDMultiply<Vector3<T1>,Vector3<T2> >{
            typedef Vector3<typename __BiOperatorType<T1,T2>::MultiType> ret_type;
            auto multiply(const Vector3<T1> &t1,const Vector3<T2> &t2)->ret_type{
                return ret_type(t1.x*t2.x,t1.y*t2.y,t1.z*t2.z);
            }
        };
        /**
         * __RTDDivide-return type deduced divide operation.
         * specialize this template to define divide operation for two customized type
         * throw std::logic_error if divided by zero
         * @tparam T1 first operand type of divide operation
         * @tparam T2 second operand type of divide operation
         */
        template<typename T1,typename T2>
        struct __RTDDivide{
            typedef typename __BiOperatorType<T1,T2>::DivideType ret_type;
            auto divide(const T1 &t1,const T2 &t2)->ret_type{
                if(t2==0)throw std::logic_error("error:cannot be divided by zero!");
                return t1/t2;
            }
        };
        template<typename T1,typename T2>
        struct __RTDDivide<Vector3<T1>,T2>{
            typedef Vector3<typename __BiOperatorType<T1,T2>::DivideType> ret_type;
            auto divide(const Vector3<T1> &t1,const T2 &t2)->ret_type{
                if(t2==0)throw std::logic_error("error:cannot be divided by zero!");
                return ret_type(t1.x/t2,t1.y/t2,t1.z/t2);
            }
        };
        template<typename T1,typename T2>
        struct __RTDDivide<T1,Vector3<T2> >{
            typedef Vector3<typename __BiOperatorType<T1,T2>::DivideType> ret_type;
            auto divide(const T1 &t1,const Vector3<T2> &t2)->ret_type{
                if(t2.x==0||t2.y==0||t2.z==0)throw std::logic_error("error:cannot be divided by zero!");
                return ret_type(t1/t2.x,t1/t2.y,t1/t2.z);
            }
        };
        template<typename T1,typename T2>
        struct __RTDDivide<Vector3<T1>,Vector3<T2> >{
            typedef Vector3<typename __BiOperatorType<T1,T2>::DivideType> ret_type;
            auto divide(const Vector3<T1> &t1,const Vector3<T2> &t2)->ret_type{
                if(t2.x==0||t2.y==0||t2.z==0)throw std::logic_error("error:cannot be divided by zero!");
                return ret_type(t1.x/t2.x,t1.y/t2.y,t1.z/t2.z);
            }
        };
        //friend overload functions of Vector3,applied to the following circumstances only:
        //Vector3(+,-,*,/)scalar, scalar(+,-,*,/)Vector3 and Vector3(+,-,*,/)Vector3
        template<typename T1,typename T2>
        auto operator+(const T1 &t1,const Vector3<T2> &t2)->typename __RTDAdd<T1,Vector3<T2> >::ret_type{
            return __RTDAdd<T1,Vector3<T2> >().add(t1,t2);
        }
        template<typename T1,typename T2>
        auto operator-(const T1 &t1,const Vector3<T2> &t2)->typename __RTDMinus<T1,Vector3<T2> >::ret_type{
            return __RTDMinus<T1,Vector3<T2> >().minus(t1,t2);
        }
        template<typename T1,typename T2>
        auto operator*(const T1 &t1,const Vector3<T2> &t2)->typename __RTDMultiply<T1,Vector3<T2> >::ret_type{
            return __RTDMultiply<T1,Vector3<T2> >().multiply(t1,t2);
        }
        template<typename T1,typename T2>
        auto operator/(const T1 &t1,const Vector3<T2> &t2)->typename __RTDDivide<T1,Vector3<T2> >::ret_type{
            return __RTDDivide<T1,Vector3<T2> >().divide(t1,t2);
        }
        template<typename T1>
        std::ostream &operator<<(std::ostream &out, const Vector3<T1> &v) {
            out<<v.x<<","<<v.y<<","<<v.z;
            return out;
        }
    }
	//for old compatibility,using class geometry::Vector3<T> in current namespace
    template<typename T>
    using Vector3=geometry::Vector3<T>;
    /**\brief encapsulate color related functions
     * some methods depend on geometry::Vector3
    */
    namespace colors
    {
        /**\brief convert rgb to hsv
          *\param[in] rgb:0~255
          *\param[out] hsv:h[0,359);s[0,1];v:0~255
          */
        void rgb2hsv(const std::vector<int> &rgb, std::vector<float> &hsv);
        /**\brief convert hsv to rgb
        *\param[in] hsv:h[0,359);s[0,1];v:0~255
        *\param[out] rgb:0~255
        */
        void hsv2rgb(const std::vector<float> &hsv, std::vector<int> &rgb);
        /**\brief note:use ColorGenerator to generate random non-repetive color instead
         * generate random color in range [pit-range,pit+range) and take the opposite of it
         * meaning generate a random color outside the given range
         * this can generate a color that contrasts to given color `pit`
         * @param pit the pit center
         * @param range the pit range
         * @return
         */
        std::vector<int> rand_color(std::vector<int> pit={0,0,0}, std::vector<int> range={64,64,64});
        /**\brief random color generator
         * used to generate random colors and assure that theuy are different
         */
        class ColorGenerator{
        private:
            std::list<geometry::Vector3<int> > history;
            int deny_range;
            int max_size;
            bool exists(geometry::Vector3<int> color);
            void rememberMe(geometry::Vector3<int> color);
        public:
            /**\brief create a random color generator
             * @param deny_range if r,g,b value of color1 are all arround(within range of
             *                   -deny_range~deny_range) r,g,b values of color2,color1==color2
             * @param max_history maximum history of appeared color(to avoid repetition);if -1=>no limit
             */
            ColorGenerator(int deny_range=5,int max_history=-1);
            /**\brief create a random color generator
             * @param background background color(this is put to history automatically)
             * @param deny_range if r,g,b value of color1 are all arround(within range of
             *                   -deny_range~deny_range) r,g,b values of color2,color1==color2
             * @param max_history maximum history of appeared color(to avoid repetition);if -1=>no limit
             */
            ColorGenerator(geometry::Vector3<int> background,int deny_range=5,int max_history=-1);
            /**\brief generate a random non-repetitive color
             * @param never_appear_again whether the generated color should never appear again
             * @return the generated random non-repetitive color
             */
            geometry::Vector3<int> operator()(bool never_appear_again=true);
        };
        //color constants
        const geometry::Vector3<unsigned char> black=geometry::Vector3<unsigned char>({0,0,0});
        const geometry::Vector3<unsigned char> red=geometry::Vector3<unsigned char>({255,0,0});
        const geometry::Vector3<unsigned char> green=geometry::Vector3<unsigned char>({0,255,0});
        const geometry::Vector3<unsigned char> blue=geometry::Vector3<unsigned char>({0,0,255});
        const geometry::Vector3<unsigned char> yellow=geometry::Vector3<unsigned char>({255,255,0});
        const geometry::Vector3<unsigned char> magenta=geometry::Vector3<unsigned char>({255,0,255});
        const geometry::Vector3<unsigned char> cyan=geometry::Vector3<unsigned char>({0,255,255});
        const geometry::Vector3<unsigned char> white=geometry::Vector3<unsigned char>({255,255,255});
        const geometry::Vector3<unsigned char> purple=geometry::Vector3<unsigned char>({128,0,128});
        const geometry::Vector3<unsigned char> orange=geometry::Vector3<unsigned char>({255,165,0});
        const geometry::Vector3<unsigned char> gray=geometry::Vector3<unsigned char>({128,128,128});
    }
    /**\brief define abbreviations for above namespaces,usage:
     * \using namespace lyxutils::abbreviations;
     */
    namespace abbreviations{
        namespace dbg=lyxutils::debug;
        namespace str=lyxutils::str_utils;
        namespace xio=lyxutils::io;
        namespace np=lyxutils::numpy;
        namespace xos=lyxutils::os;
        namespace geo=lyxutils::geometry;
        namespace xcolor=lyxutils::colors;
        //Vector of different common types
        using Vec3u=lyxutils::geometry::Vector3<unsigned char>;
        using Vec3i=lyxutils::geometry::Vector3<int>;
        using Vec3f=lyxutils::geometry::Vector3<float>;
        using Vec3d=lyxutils::geometry::Vector3<double>;
    }
}
#endif // !UTILS_H